%% 计算单元能量（动能ke(i)和势能se(i)）
% 某一时刻，全部单元。

function [se,ke]=func_energy(m,k,nBall, a, da)

% 单元势能（两个“半Spring”的势能）
se=zeros(nBall,1);
se(1)=0.5*0.5*(2*k)*(a(2)-a(1))^2;
se(1)=se(1)+0.5*0.5*k*(a(2)-a(3))^2;
for i=2:1:nBall-1
    se(i)=0.5*0.5*k*(a(i+1)-a(i))^2;
    se(i)=se(i)+0.5*0.5*k*(a(i+1)-a(i+2))^2;
end
se(nBall)=0.5*0.5*k*(a(nBall+1)-a(nBall))^2;
se(nBall)=se(nBall)+0;

% 单元动能（Node的动能）
ke=zeros(nBall,1);
for i=1:1:nBall
    ke(i)=0.5*m*da(i+1)^2;
end

end 